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Abstract— The paper presents an application of adaptive finite dement methods to the modeling of 
low-cycle continuum damage and life prediction of high-temperature components. The major objective 
is to provide automated and accurate modeling of damaged zones through adaptive mesh refinement and 
adaptive time-stepping methods. The damage modeling methodology is implemented in an usual way 
by embedding damage evolution in the transient nonlinear solution of elasto-viscoplastic deformation 
problems. This nonlinear boundary-value problem is discretized by adaptive finite element methods. 
The automated h -adaptive mesh refinements are driven by error indicators, based on selected principal 
variables in the problem (stresses, nonelastic strains, damage, etc.). In the time domain, adaptive time- 
stepping is used, combined with a predictor-corrector time marching algorithm. The time step selection 
is controlled by required time accuracy. In order to take into account strong temperature dependency 
of material parameters, the nonlinear structural solution is coupled with thermal analyses (one-way 
coupling). Several test examples illustrate the importance and benefits of adaptive m es h refinements in 
accurate prediction of damage levels and failure time. 


L INTRODUCTION 

Unified viscoplastic constitutive models have evolved 
over the last 20 years to provide a means for 
analytically representing a material response from the 
elastic through the plastic range, including strain-rate 
dependent plastic flow, creep and stress relaxation. 
These theories are guided by physical considerations 
including dislocation dynamics and are based on the 
principles of continuum mechanics. 

The first multi-dimensional formulations of elasto- 
viscoplastic constitutive equations was due to Bodner 
and Partom [1]. Since then, a number of constitutive 
models have appeared; many of these theories are 
summarized in review articles that appear in Ref. [2]. 
A NASA-Lewis Research Center sponsored research 
program (HOST) [3, 4] conducted by the Southwest 
Research Institute recently concluded a four year 
research effort to develop unified constitutive models 
for isotropic materials and to demonstrate their use- 
fulness for analysis of high temperature gas turbine 
engines. One result of this study is material property 
data for high temperature nickel-based alloys over 
a wide temperature range. The unified models 
employed were those of Bodner-Partom and Walker. 
More recently, a viscoplastic constitutive model was 
developed at NASA Lewis by Freed [5] and Freed 
and Verilli [6], and then refined by Freed and Walker 
[7]. Relevant material constants were determined 
experimentally for aluminum, copper and nickel. 


Unified viscoplastic theories have been implemented 
by a number of finite element researchers. The first 
use of unified constitutive modeb in a finite element 
program was by Newman et aL\$l Later, under 
the NASA HOST program, the Bodner-Partom and 
Walker modeb were implemented into finite element 
codes by General Electric and Pratt and Whitney. 
Kaufman et ai. [9] describes these efforts as applied 
to gas turbine components. Moreno and Jordan [10], 
under the NASA HOST program, developed and 
applied a unified constitutive model to gas turbine 
combustors using the MARC code. Recently, a 
finite element viscoplastic study of cylindrical thrust 
chambers was presented by Arya [11]. Detailed studies 
of several rate-dependent plasticity modeb and their 
numerical implementation using adaptive finite ele- 
ment methods were conducted by Bass and Oden [12]. 
These studies were extended for thennoviscoplastic 
analysis by Thornton et ai. [13]. 

Modeb of material damage and reliability studies 
are relatively recent developments, but significant pro- 
gress has been made in understanding and modeling 
these problems. The modeb of continuum damage 
development, which are of primary interest here, 
are based on micro-mechanical studies of micro-crack 
nucleation and growth, discussed by Ashby and Raj 
[14], Greenwood [15] and many other researchers. 
Phenomenological descriptions of these models are 
usually based on the early work of Kachanov [16], 
where the general evolution law was introduced for 
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Fig. 1. Adaptive mesh refinements for damage modeling, 
(a) resolution of distributed damage zones, (b) localization 
of macro-crack initiation. 

damage parameter. These, models fit very well into 
the framework of unified viscoplastic theories, and 
several detailed evolution forms have been proposed 
(see, for example, papers of Bodner and Chan [17], 
Leckie and Hayhurst[18], Simo[19], and references 
therein). The equations presented in these works are 
convenient for the description of low-cycle damage, 
where the damage occurs relatively quickly under 
a general stress history. For problems with cyclic 
loads, the evolution of damage is more conveniently 
expressed in terms of the number of cycles. Theories 
of this type were reviewed by Nelson [20]. In high- 
cycle fatigue modeling, they are computationally 
more efficient than time integration methods. To 
date, no adaptive finite element solutions have been 
reported for continuum damage modeling in the 
literature. 

Therefore, the main objective of this study was 
to prove the feasibility and practical benefits of using 
adaptive finite element methods to model continuum 
damage and to improve life prediction capabilities for 
high-temperature components and other structures 
subjected to hostile thermo-mechanical environments. 

The major benefit of such an application of adaptive 
methods is that the adaptive mesh refinement based 
on rigorous error estimators, can provide automated 
and exact resolution of critical damage zones (Fig. la) 
and, moreover, an automated indication of possible 
indication of macro-cracks (Fig. lb). Independent of 
adaptive mesh refinement, an adaptive time-stepping 
method can provide automated control of the accur- 
acy of time-marching procedures, and greatly minimise 
the computational effort necessary to obtain a final 
solution. 

2. UNIFIED VISCOPLASme THEORIES AND 
CONTINUUM DAMAGE MODELINGS 

Complex mechanical behavior of many high- 
temperature engineering materials can be improved 
by the constitutive models known as unified elasto- 
viscoplastic theories. These theories, consistent with 
micro-mechanics of metallic materials, describe in a 
unified manner various aspects of elastic and non- 
elastic behavior of structures. More, these approaches 
are readily expandable to include such effects as crack 


nucleation, continuum damage, etc. In this study the 
Bodner-Partom [1, 17] and the Freed-Walker con- 
stitutive theories [5-7] were implemented for damage 
modeling. It should be noted that presently there exists 
a wide variety of viscoplastic constitutive theories 
and corresponding damage models. However, the 
primary objective of this paper is not to study these 
models, but to develop and demonstrate adaptive 
finite element techniques in application to damage 
modeling. Hence the above two constitutive theories 
were implemented in this study. Importantly, the 
adaptive finite element techniques developed herein 
is applicable to a much wider class of continuum 
damage models. 

2.1. Basic formulations 

The intemal-state-variable constitutive formulations 
are to be used in the solution of infinitesimal, quasi- 
static, isothermal, initial-boundary-value problems. 
Consider a viscoplastic structure occupying a region 
Q with boundary dQ. The behavior of the structure 
is described by the following system of differential 
equations: 


(1) Equilibrium in rate form: 

* w + 6,«0, (1) 

where denote components of the stress 
tensor, b t are the body force components per 
unit volume, and the summation convention is 
employed. 

(2) Kinematic relation for velocity gradients: 

+ + ( 2 ) 


where denotes the total strain components 
and superscripts E and P denote elastic and 
inelastic strain components, respectively. The 
components of the displacement rates are ii,. 

(3) Constitutive relations: 


f 5 -/#(**. Z k ) (no sum) 
it - */(**» Z k \ (no sum) J 


( 3 ) 


where represent Hooke’s tensor of elasticity 
parameters, 0 *, are components of a tensor of 
thermal expansion parameters, and t represents 
the rate of the change in temperature. Both 
and a*/ are temperature dependent. The consti- 
tutive functions are f ¥ and g t where Z k represent 
internal state variables. These functions and 
state variables characterize the viscoplastic 
response of the material. 


The description of the problem is completed by 
prescribing the boundary and initial conditions. 
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(4) 


where dQ^dfl, uq&j, 'ft, are prescribed surface 
displacement rates, the components of a unit 
normal vector, and ft, the prescribed surface traction 
rates. The initial conditions indude specifying the dis- 
placements, stresses, and internal state variables, i.e. 


w, (*, 0), ofa OX Z,(x, 0) xeQ (5) 


2.2. Continuum damage models 

The continuum damage modeling describes the 
deterioration of mechanical properties of materials 
by development of intergrain cracks and voids. 
According to material science theories presented, for 
example, by Ashby and Raj [14] or Greenwood [15], 
the growth of damage at grain boundaries is expressed 
in terms of two mechanisms: nucleation and growth. 
The damage parameter to is then introduced, and can 
be interpreted as the ratio of damaged cross section 
(-4 vow) (or, alternatively, volume) to the original area 
): 



The micro-cracks considered here are in the range 
of 0.01 mm in length. The rupture criterion is to — 1, 
which corresponds to the saturation of the material 
with voids. Alternatively, a single crack may grow 
to a size on the order of 1 mm. In the latter case, 
the crack is too big to be treated in a continuum 
sense, and its propagation should be followed using 
the methods of fracture mechanics. 

Continuum damage modeling is easily incorporated 
into the viscoplastic constitutive theories by intro- 
ducing the evolution law for the damage parameter. 
Usually the general Kachanov form based on 
separation of variables is accepted: 

d> (7) 

where a is the stress state in the body. 

More generally, other parameters can be used in 
the evolution law. Moreover, the damage parameter 
is often decomposed into isotropic and anisotropic 
components. Within this class, two low-cycle damage 
models have been considered in this study based on 
the Bodner and Chan [17] and Greenwood [15] models. 
However, there are several problems with practical 
application of directional damage, reliability of these 
models and conducting experiments relevant for the 
evolution of necessary parameters. Even the extensive 
experiments presented in Refs. [3,4,17] did not 
provide all the necessary data and, hence, the damage 
model was limited to the isotropic damage. There- 
fore, numerical implementation of this study was also 
limited to isotropic damage. 


Once the damage parameter o> is calculated, it 
should be somehow incorporated into the constitutive 
theory. This is based on a straight forward interpret- 
ation of the damage parameter, co, as a ratio of area 
of voids in a given cross-section to the original area, 
as shown in eqn (6). The result of this reduction of 
area is that the real stress in a material is higher 
th an the nominal stress 

4* ml 



In the Bodner-Partom theory this results in the 
modification of the nonelastic strain rate equation. 
The details of this modification can be found inf Ref. 
[17]. The same principle, applied to Freed-Walker 
theory [7], leads to the use of a real stress rather than 
a nominal stress in the evolution equation. Therefore, 
to incorporate continuum damage, the stress deviator 
(Sy) in the viscoplastic model used in this study is 
being replaced by: 

S ¥ -Z^(\-w)S ¥ . (9) 

This replacement means that the stresses stored in 
the program are nominal stresses, but the nonelastic 
flow rule operates on real stresses. Consequently, the 
back stress B and the drag stress D in their model are 
real stresses. 


X FINITE ELEMENT SOLUTION OF VTSCOPLASTIC 
PROBLEMS WITH DAMAGE 


The boundary-value problem described in the 
previous sections is strongly nonlinear and can only 
be solved by numerical techniques, such as the finite 
element method. Application of this method requires 
recasting of the problem in a weak, variational form. 


3.1. Weak formulation 

The rate form of the equilibrium equations de- 
scribed earlier can be viewed as a constraint condition 
on the constitutive equations. A weak formulation 
of the equilibrium equations, sufficient to include 
the necessary constraint condition, is developed here. 
The equilibrium conditions in terms of the strain rates 
are of the form 


where 


► 

(no sum) 


( 10 ) 
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where a y and Z, are understood to depend upon u Q 
and to satisfy initial conditions. 

In order to obtain a weak formulation of a 
boundary-value problem the space of displacement 
rates (solution space of trial functions) is introduced. 


U-{ue[W m ''(Q)Y', i(x) — ► 0 as || x || — ► ao} (12) 


Conversely, any solution of the governing equations 
and boundary conditions will also satisfy eqn (15). 

3.2. Finite element formulation 

A finite element approximation of eqn (15) is 
obtained by introducing approximations of the trial 
and test functions: 


**(*)- 2 »/*) 

. . x-"'" 

«»(*. 0 - £ «y(< w*). 


ciA 


(16) 


and the space of test functions 

F«{v6 ▼(*)—► 0 as ||x J — ► co}, (13) 

where O is a computational domain, N is the 
dimension of the physical space (2 or 3), and 
is the Sobolev space, where specific values of m, p 
and q depend on the particular forms of constitutive 
equations governing the material under consideration. 

Multiplying the equilibrium equation by a suitably 
smooth test function r, — ^(x) and integrating over 
Q, the weak form of the rate equilibrium equations 
can be obtained: 


where N is the number of nodes, i indicates the vector 
components, is a basis function, and uj, Jj and «J(0 
are die values of u[ t and u{(t) node j at time r, 
respectively. 

Recasting these approximations into a matrix form 
and employing standard notation, gives the time 
element approximation for u(jc, /) at time / 

| B T EBa dfl = J B T K f dfl + j* B T Eaf dn 

+ j * T fcdfl+ [ (17) 

Jo J**, 





« 0 for all r, in K. 


(14) 


After the substitution of the constitutive equations, 
application of the divergence theorem and a grouping 
of terms the following variational problem is obtained: 
Find a displacement rate field t — ► u(x, /) in V + {«} 
such that for every t 


j* EyuUtjVy dQ » J* Ej£,lu r <jdO 

+ J 

+ j* &iV t d Cl + J i t Vi ds 

for all v t in V. 


(15) 





The rates of nonelastic strains i f u are obtained from 
the relevant constitutive theory, and the pointwise 
values of rates of temperature t were read in from 
the file created by thermal analysis program. Here 
dn and dr are volume and surface measures, u is 
any function defined on Cl such that its trace on the 
boundary segment 30, isk (where u,3Q, =*-£,% It is 
easily verified that any sufficiently smooth solution of 
eqn (1) will also satisfy the governing equations. 


where O denotes the element volume, dOj denotes 
an element surface where tractions are defined, and 
B is the strain displacement matrix. 

The temperature affects the viscoplastic structural 
analysis directly in three ways: (1) the elasticity 
matrix E(r) and the coefficient of thermal expansion 
a(T) depend on temperature, (2) nodal toads-f^T") 
depend on the local temperature rates and (3) 
several parameters in the viscoplastic constitutive 
model are temperature dependent 

3.3. Viscoplastic solution method with damage modeling 

In this section, an algorithm used for the solution 
of viscoplastic problems with damage modeling is 
presented. In addition to this basic algorithm, an 
interface was provided for the thermal analysis, 
wherein temperature histories can be read into the 
structural code. 

An uncoupled (sometimes called one-way coupled) 
formulatoin was used for this approach, wherein 
the thermal problem is solved first followed by the 
viscoplastic analysis. The transient thermal problem 
is solved by time marching with a time step A t T , and 
the nodal temperatures at successive times *, , t 2 , . . . 
are obtained. These temperature vectors arc used as 
input to the structural analysis. 

The first computation in the structural analysis is 
to solve an initial statics problem if the initial tem- 
perature distribution T(x, 0) is not equal to a uniform 
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reference temperature or if any initial static loads are 
present. The results of this analysis arc the initial con- 
ditions (displacements and stresses) for the transient 
viscoplastic analysis. 

The viscoplastic analysis time marches with a 
time step A /, . Experience has shown that the time step 
required for the structural analysis is usually smaller 
than for the thermal analysis, i.e. A T % < A/ T . At inter- 
mediate times in the structural analysis, the tempera- 
tures are linearly interpolated from the temperatures 
known at the beginning and end of the larger thermal 
time intervals. 

A computationally efficient method lor tracing the 
solution of inidal-boundary-value p robl ems in visco- 
plasticity, which utilize a hidden variable constitutive 
formulation is described here. The strategy employed 
in the viscoplastic algorithm is as follows: 

With the initial distribution of stress, temperature, 
and internal variable specified, use the equilibrium 
equations to obtain the nodal displacement rates. 
Then integrate the constitutive equations forward in 
time at the element Gauss integration points. With 
updated values of the stress, temperature, and internal 
variables at the new time, the equilibrium equations 
are solved again. This sequence is determining the 
nodal displacement rates, then advancing the consti- 
tutive equations in time is continued until the desired 
history of the mitial-boundary-vaiue problem has 
been obtained. Thus, the algorithm proceeds through 
the following steps: 

(1) At time /„ * 0, initialize a 9 , Z, and m, for each 
element. 

(2) Calculate -/*(**, Z*,m,) Cor each element 
at t — [thus determining the right-hand side 
of eqn (17)]. 

(3) Assemble and solve the equilibrium condition 
for u*(fj. 

(4) Calculate i 9 (t M ) - Bu*(/J for each element. 

(5) Calculate &Jt 9 ) - !#,(/.) -fJ(O) - 

(6) Calculate £(t m ) - gt(p 9 (t m \ (/.)). 

(7) Calculate 

(8) Integrate % and d> { forward over some 
appropriate A/, for each dement to get 

(9) If 1 - t 9 \ At , < t ** goTo 2, otherwise st5p. 

Step 8 characterizes an explicit scheme, and this 
step can be replaced by a subroutine for implicit 
method if needed. If predictor-corrector type inte- 
gration schemes are used, a slight modification is 
necessary beginning with Step 8 to be described in 
Section 4.2. The computational method above has 
been presented for a constant time step At,. The com- 
putational experience by several investigators [9, 12, 
16, 21] indicates that a very smal time step can be 
required because of the mathematically stiff nature of 
the constitutive equations describing the internal state 
variables. To gain improved efficiency and reliability, 
a variable time step algorithm has been implemented 


in this study. This technique will be described in the 
next section. 

4. ADAPTIVE COMPUTATIONAL TECHNIQUES 

In this study, adaptive computational techniques 
were used to provide the maximum reliability of 
results at the minimum computational cost. These 
techniques are used to adapt their parameters to 
changing characteristic of the solution during the 
solution process. In this study, two bask adaptive 
methods were used, namely h -adaptive mesh refine- 
ment and adaptive time stepping. The basics of 
these techniques will be discussed on the following 
subsections. 

4.1. Error estimation and adaptive mesh refinement 

Adaptive mesh refinement, based on rigorous 
error estimates, enables to resolve a bask question 
of modern computational mechanics- what is the 
accuracy of the numerical solution? equivalently, 
how good are the answers? To systematically reduce 
error by adaptively changing mesh sizes or spectral 
orders of approximations, one must obviously have 
some means to judge the distribution of error in a 
numerical solution. 

Some of the most important results related to error 
estimation and adaptive techniques for finite element 
methods and details of their derivations and proofs 
can be found in Refs. [22-24]. The error estimators 
and indicators are the basis for adaptive mesh refine- 
ments, which modify the mesh in order to minimize 
the error at minimal computational cost. In this 
study, interpolation error estimators combined with 
h -refinement techniques have been implemented for the 
solution of dasto-viaooplastk problems with damage. 
The interpolation methods use the interpolation theory 
of finite elements in Sobolev norms to produce rapid 
(sometimes crude) estimates of the local error over 
individual elements. Simplifying this somewhat gives 
an error estimate HefineH as the maximum difference 
of the effective strain within the union of that par- 
ticular element Q, and its neighboring dements. The 
h -methods are the mesh refinement methods, based 
on local subdivision of dements into smaller elements. 
An h-type refinement strategy may be outlined as 
follows [11]: ^ 

(1) Given a solution at time r*, solve the initial- 
boundary-value problem until time r“ +l and calculate 
the estimated error for each element. The interval 
/"♦ 1 2 3 4 5 6 7 8 9 — /• is much larger than the time inarching step 
Ar. 

(2) Determine the magnitude of the largest elemental 
error and set a threshold value for refinement 
a£^u, where a is a fraction between zero and one. 

(3) Determine the number of demeats to be refined 
using the threshold value. If this number of elements 
is greater than fiN 9 where N t is the number of ele- 
ments and fi is a fraction between zero and one, a 
uniform distribution of error exists and no refinement 
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is made (go to Step 1). Otherwise, refine all elements 
for which E k >aE imx9 where ^ is the error for a 
typical element A: in the mesh. 

(4) If refinement has occurred, go back to time 
interpolate the solution on a new mesh, and resolve 
the boundary value problem on the new mesh. Then 
again estimate the elemental error at time /" +1 . 

(5) If element errors exist such that E k > ot E^ go 
to Step 3. Otherwise the errors have been reduced 
sufficiently, go to Step 1 with a new refined mesh, and 
new starting time /* +l . 

Ifis quite essential for the accuracy of the solution 
that after mesh refinement at time the solution 
process is backed to the time t ■ This assured that the 
interval is covered on adequately refined 

mesh, so there is no excessive accumulation of error. 

4.2. Adaptive time stepping techniques 

Time-dependent problems are usually solved by 
time marching techniques, where the time domain is 
discretized into a number of steps f, , t 2t . . . , t„ and 
the solution proceeds by solution of a sequence of 
incremental problems. In the simplest algorithms, 
the time step interval is kept constant throughout 
the whole computational process. In most practical 
applications, a fixed time step leads to relatively 
inefficient solutions, and adaptive adjustments of At 
are very desirable. In this section, an adaptive time- 
stepping technique for thermo-viscoplastic structural 
analysis with damages is described. 

For simplicity, consider the single ordinary differ- 
ential equation, 

y-f(y*t). (18) 



Fig. 2. Adaptive timestepping algorithm for viscoelasto- 
plastic evolution problems. 


variables o) h with the maximum of these components 
selected as the controlling error. The corrector phase 
is the modified Newton scheme, 



(22) 

y?+i,=y, 

(23) 


A flow chart depicting the adaptive scheme is 
shown in Fig. 2. 


The solution is advanced using a predictor- 
corrector scheme. The predictor phase consists of an 
Euler step: 


yt+A,-y,+Aty, 09 ) 

+ A/). (20) 

An error indicator £[12,22] is then computed 
from 


£ 


vr+j 


( 21 ) 


The error indicator is next compared with a preset 
error criterion and if the criterion is met, the timestep 
is small enough to proceed to the corrector stage. 
Otherwise, the predictor phase for eqs (18) — (21) is 
repeated with a smaller time step. On the other hand, 
if the error is smaller than the user-defined value E mimf 
the time step is increased to improve computational 
efficiency. For the viscoplastic evolution equations 
with damage, the control variables used to calculate 
the error indicator were the components of a stress 
tensor <fy, internal state variables Z„ and the damage 


5. NUMERICAL EXAMPLES 

In this section, selected numerical examples oriented 
on verification of damage modeling and on illustration 
of advantages of adaptive methods in application 
to damage modeling are presented. In particular, 
discussed are: (1) a simple tension creep test with 
experimental verification; (2) an isothermal structural 
problem with adaptive mesh refinement; (3) a thermo- 
structural problem corresponding to the space shuttle 
main engine (SSME) combustion chamber liner. 

5.1. Tensile creep test 

The main purpose of this example is to verify the 
viscoplastic constitutive model and damage modeling 
by comparison of numerical predictions with avail- 
able experimental data [17]. The example is a tensile 
creep test of nickel superalloy B1900 + Hf subjected 
to constant stress level of 215 MPa and temperature 
of 982°C. For this case, both tertiary creep corre- 
sponding to continuum damage and final rupture 
of the specimen were observed. The resulting time- 
strain curves and damage histories are 'presented 
and compared with experimental data in Fig. 3. For 
comparison, the same examples were solved without 
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damage modeling. The results of these comparisons 
are presented in Fig. 4. 

The following conclusions can be drawn from 
this example and other tests shown in Ref. [17]: 

(1) the damage evolution models describe correctly 
the tertiary creep zone and transition to rupture; 

(2) the evolution of damage and rupture time are 
extremely sensitive to the stress level — small variation 
of stress or material constants for damage cause 
considerable change in the rupture time. A more 
extensive study of this issue is presented in Ref. [17]; 

(3) the actual rupture of the experimental specimen 
occurs at damage levels smaller than the theoretical 
rupture condition (co = 1). This is most probably 
caused by the fact that local damage concentrations, 
say on material inhomogeneities, cause localized 


growth of continuum damage, transition to macro- 
cracks, and their propagation before a uniform satur- 
ation of the whole cross-section with micro-cracks 
can occur; (4) without damage modeling, the tertiary 
creep zone and rupture of the material are not 
predicted by the viscoplastic constitutive theories — 
thus damage modeling is necessary to capture these 
phenomena. 

5.2. L -shaped domain 

A plane stress L -shaped domain subjected to pre- 
scribed displacement rates is presented as an example 
of a mesh-adaptive isothermal analysis with damage 
modeling (Fig. 5). In particular, displacements on 
boundary segments indicated in Fig. 5 grow during 
10 s to a value of —0.001, and then are held steady 


/ Material: B1900 + Hf 
/ / Temperature: 1000° C 


Problem statement: 



u* = -0.001 




Load history: 



CAS— MS 059 




until the final time of 1000 s. The temperature is equal 
1000°C, and the material is a superalloy B1900 + Hf. 
The problem was solved with adaptive mesh refine- 
ment. A comparison of initial mesh and the final 
adapted mesh is shown in Fig. 6. Contours of stress 
intensity and of damage are presented in Fig. 7. As 
expected, the adaptive procedure correctly refined the 
mesh around the comer singularities. The adapted 
mesh identified the maximum damage zone and 
predicted a high damage level, Q)«». * 0.45. In fact, 
the damage zone is so concentrated that initiation 
of macro-crack can be predicted here. Although the 
location of crack initiation (comer of the L -shaped 
domain) is obvious in this test problem, in more 
general practical situations the appropriate adaptive 
procedure can also correctly identify such cases. 

5.3. SSME combustion chamber liner 

An adaptive analysis of an element of the SSME 
combustion chamber liner, subjected to thermal and 
mechanical loads has been taken as an illustration 
of the application of damage modeling to practical 
problems in aerospace structures. The problem for 
this example is illustrated in Fig. 8. The cross section 
of one cell of the linear is almost rectangular, with the 
internal coolant passage subjected to a pressure of 


Fig. 7. Effective stress contourr (a) and damage level (b) 
for an L-shaped domain-refined mesh. 


3900 psi (25.5 MPa). The initial temperature is of the 
order of — 200°C, with almost uniform distribution 
throughout the cross-section. The thermal load history 
corresponds to a rapid heating of the liner from the 
inside in 5 s up to a temperature close to 600°C. 
This tem p erature gradient is held almost constant 
for 10 s, and then the liner is cooled down in 5 s. 
An idealization of this history is shown in Fig. 8. 

The above problem was solved with boundary 
conditions, corresponding to a single cell of an axi- 
symmetric liner, as shown in Fig. 9. The initial mesh 
for this problem is shown in Fig. 10. Starting from 
this mesh, a finite element adaptive procedure gener- 
ated the refined mesh, as shown in Fig. 11. As will be 
shown later, this mesh correctly identifies regions of 
high stress and damage. For this problem, a combin- 
ation of internal pressure and temperature gradient 
causes an extensive deformation of a bottom wall 
shown with exaggerated scale in Fig. 11. The corre- 
sponding temperature distribution is shown in Fig. 12 
— high gradients are present near the bottom wall of 
the cell. The high temperature gradient and internal 
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Geometry: 



0.040 in 

Material: copper 
Thermal history: 




Fig. 8. SSME combustion chamber liner: problem statement 
and simplified thermal load history. 



Fig. 9. Boundary conditions for the combustion chamber 
liner. 


pressure cause high stress levels and some plastic 
deformation of the bottom wall as shown in Figs. 13 
and 14. Due to the mostly compressive nature of 
stresses at this stage, damage development is very 
low and does not exceed co « 10~ 7 at the end of the 
heating period. 

During a cooling phase the situation is generally 
reversed — due to the specific boundary conditions 
and the plastic yield which occurred during heating, 
high tensile stresses develop in the bottom wall of the 
cell. The deformed configuration at the end of the 
cooling period is shown in Fig. 15. One can observe 
a thinning of the bottom wall caused by a combin- 
ation of tensile stresses and internal pressure. The 
distribution of stresses at the end of the cooling 
period is shown in Fig. 16 and the corresponding 
principal plastic strain is illustrated in Fig. 17. In the 
cooling phase, the damage development is accelerated 



Fig. 10. SSME liner-initial mesh. 



Fig. 11. SSME liner-deformed configuration at the end of 
the heating period: / — 15 seconds. 


in the bottom wall — the corresponding distribution 
of damage is shown in Fig. 18. 

To better illustrate the development of the damage 
pattern, histories of damage for the inner lower 
corner and for the center of lower face of the cell are 
shown in Fig. 19. It can be 9 een that damage develop- 
ment is indeed accelerated during the cooling phase. 
One can expect that further damage will evolve in the 
consecutive cycles of heating-cooling. This example 
clearly illustrates the potential of adaptive finite 
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TIME (SEC) X 10—1 


a) 


rsi 



b) 


Fig. 19. SSME liner-histories of damage at various locations in a cross-section. Note a different scale for 

each plot. 


element methods in the identification of critical zones 
and damage patterns of aerospace structures. 

& CONCLUSIONS 

The formulations and results presented in this 
paper clearly confirm the feasibility of developing 
adaptive finite element methods for modeling damage 
evolution and for making life predictions of high- 
temperature structures. 


The error estimation and adaptive mesh refinement 
not only provide accuracy control and an optimal 
solution at a minimum cost, but are also capable of 
pinpointing locations, where concentrated continuum 
damage growth leads to initiation of macro-cracks. 

Similarly, the adaptive time-stepping guarantees 
error control and optimal adjustments of time steps, 
with the actual values of the time steps often varying 
by the factors of 1000 depending on the evolution of 
the nonlinear deformation in the structure. 
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As far as damage modeling is concerned, the 
following remarks are applicable: (l) while theoretical 
formulations of the damage modeling are presently 
quite mature, the number of reliable experimental 
studies dedicated to the determination of material 
constants and the verification of theory is still very 
limited. (2) The rate of continuum damage develop- 
ment is very sensitive to the temperature and the 
stress level — a small variation of these parameters 
leads to considerable changes in the predicted time to 
rupture. This sensitivity emphasizes the importance 
of an accurate prediction of the stresses by the 
application of error estimation and mesh refinement. 
(3) This paper addressed only the low-cycle damage 
development, wherein rupture happens after a 
relatively short time, and the damage parameter can 
be effectively integrated as a function of time. (4) For 
many practical applications it is also important to 
consider high-cycle damage, wherein it is much more 
efficient to integrate damage over the number of cycles, 
rather than directly over time. The extrapolation 
techniques are useful in cases of long-term damage 
evolutions. 
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